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Abstract 

Quantitative real-time polymerase chain reaction (qPCR) is an efficient and widely used 
technique to monitor gene expression. Housekeeping genes (HKGs) are often empirically 
selected as the reference genes for data normalization. However, the suitability of HKGs used as 
the reference genes has been seldom validated. Here, six HKGs were chosen {actin A3, actin Al, 
GAPDH, G3PDH, E2F, rp49) in four lepidopteran insects Bombyx mori L. (Lepidoptera: 
Bombycidae), Plutella xylostella L. (Plutellidae), Chilo suppressalis Walker (Crambidae), and 
Spodoptera exigua Hiibner (Noctuidae) to study their expression stability. The algorithms of 
geNorm, NormFinder, stability index, and ACt analysis were used to evaluate these HKGs. 
Across different developmental stages, actin Al was the most stable in P. xylostella and C. 
suppressalis, but it was the least stable in B. mori and S. exigua. Rp49 and GAPDH were the most 
stable in B. mori and S. exigua, respectively. In different tissues, GAPDH, E2F, and Rp49 were 
the most stable in B. mori, S. exigua, and C. suppressalis, respectively. The relative abundances 
of Siwi genes estimated by 2~ AACt method were tested with different HKGs as the reference gene, 
proving the importance of internal controls in qPCR data analysis. The results not only presented 
a list of suitable reference genes in four lepidopteran insects, but also proved that the expression 
stabilities of HKGs were different among evolutionarily close species. There was no single 
universal reference gene that could be used in all situations. It is indispensable to validate the 
expression of HKGs before using them as the internal control in qPCR. 
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Introduction 



Estimating transcript expression levels across 
different conditions is an important aspect of 
work in molecular biology. Many techniques 
such as northern blot, microarray, and 
quantitative real-time polymerase chain 
reaction (qPCR) have been developed to 
determine gene expression levels. Among 
these, qPCR is one of the best choices and is 
widely applied in monitoring transcript-level 
changes in gene expression. 

In several steps of the qPCR protocol, 
including preparing material, RNA 
purification, cDNA synthesis, and PCR 
procedures, non-specific variations may be 
introduced. For example, the efficiencies of 
reverse transcription (RT) and PCR may be 
varied between samples. Therefore, selection 
of one or more reference genes, also called 
internal controls, is very important to 
eliminate the inter-sample variations in 
analyzing qPCR data. An ideal reference gene 
is constantly expressed in all samples under 
certain condition. Housekeeping genes 
(HKGs) are thought to have key roles in 
cellular processes and are expressed at a 
constant level, and thus are often used as 
reference genes in qPCR. There are two 
strategies to choose a reference gene. One is 
to use HKGs as reference genes that have 
already used in previous studies; another is to 
use HKGs that are homologous to a widely 
used reference gene in other model species. 

These empirical choices often lack 
experimental support. Increasing evidence has 
suggested that HKGs are not always 
expressed stably in all experimental 
conditions. The a-actin gene, which has been 
widely used as a reference gene, was reported 
to be highly regulated by matrigel, and 
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therefore unsuitable as an internal control in 
this condition (Selvey et al. 2001). Several 
housekeeping genes, such as a-tubulin and 
GAPDH, were also reported to have varied 
expression levels in some instances such as 
different tissues or different developmental 
stages (Brunner et al. 2004; Trivedi and Arasu 
2005; Pei et al. 2007; Dhar et al. 2009; Dong 
et al. 2010; Lord et al. 2010; Wang et al. 
2010). Some studies also proved that choosing 
unsuitable reference genes produced low 
precision or misleading results (Jian et al. 
2008). Therefore, it is necessary to validate 
the expression stabilities of HKGs when 
choosing them as reference genes in qPCR. 

Though qPCR has been widely used to 
estimate the gene expression in insects, very 
few studies have been performed to validate 
the expression stability of HKGs. In Bombyx 
mori L. (Lepidoptera: Bombycidae), 
translation initiation factor 4A, translation 
initiation factor 3 subunit 4, and translation 
initiation factor 3 subunit 5 were reported to 
be the most reliable reference genes during 
metamorphosis (Wang et al. 2008). However, 
the most commonly used reference gene, 
cytoplasmic actin, varied drastically 
throughout metamorphosis development 
(Wang et al. 2008). Ribosomal protein genes 
RPS3, RPS18, and RPL13a were the most 
stable genes in Tribolium castaneum exposed 
to Beauveria bassiana (Lord et al. 2010). 
Here, four economic and agriculturally 
important lepidopteran insects were chosen, 
including B. mori, Plutella xylostella L. 
(Plutellidae), Chilo suppressalis Walker 
(Pyralididae), and Spodoptera exigua Hiibner 
(Noctuidae) to evaluate the expression 
stability of six HKGs {actin A3, actin Al, 
GAPDH, G3PDH, E2F, and rp49). Reliable 
reference genes are proposed in these four 
lepidopteran insects. 
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Materials and Methods 

Sample collection 

All insects were maintained at 25 ± 1 °C 
under a photoperiod of 16:8 L:D and 80% RH. 
The silkworms (B. mori) were fed on 
mulberry leaves, diamondback moths (P. 
xylostella) were reared on radishes, rice stem 
borers (C. suppressalis) were fed on rice 
seedling, and beet armyworms (S. exigua) 
were fed on an artificial diet (Huang et al. 
2002; Merkx-Jacques and Bede 2005). The 
different developmental stages collected were 
egg, larvae (collected at the first day of each 
instar), pupa, and adult, for each insect. 
Tissues from the head, midgut, ovary, testis, 
fat body, malpighian tube, and epidermis were 
dissected from final instar larvae and kept at 
-70 °C for RNA purification. The silk glands 
of the B. mori were also collected for 
experiments. Two sets of samples were 
collected for experimental replication. About 
30 insects were collected for each sample, and 
all experiments for each set of samples were 
repeated three times. 

Collecting the sequences of candidate 
reference genes 

A list of housekeeping genes (HKGs) was 
compiled that have been commonly used as 
the internal controls for normalization in 
qPCR data analysis in four insects by 
reference mining. The sequences of selected 
HKGs of P. xylostella, C. suppressalis, and S. 
exigua were obtained by searching 
transcriptome data. The sequences of B. mori 
genes were downloaded from the refseq 
database of NCBI. All the primers were 
designed using Beacon Designer 7 (Premier 
Biosoft, www.premierbiosoft.com). The 
sequences of primers are given in Table 1. 
Siwi-l and Siwi-2 genes were chosen to study 
the impact of reference genes on the qPCR 
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data analysis. The forward primer sequence of 
Swi-lwas 5'- 
GTGCGGGTCTGCCTGAGTTG-3 ' and the 
reverse primer sequence of Siwi-l was 5'- 
GGTGTGCTTGTGAATCTCCTGTTG-3 ' . 
The forward primer sequence of Siwi-2 was 
5 ' -AGCCTACCCTGACCTATGTTG-3 ' and 
the reverse primer sequence of Siwi-2 was 5'- 
ACCAGTCCCGTCTCGTTATAC-3 ' . 

RNA purification 

The samples were frozen with liquid nitrogen 
and homogenized in a tissue grinder. Then, 1 
mL TRIzol reagent (GIBCO, 
www.invitrogen.com) was added to 
homogenized samples. Total RNA was 
isolated following the recommended 
procedures. Genomic DNA was removed 
from total RNAs by DNase I treatments 
following the protocol of the DNA-free kit 
(Ambion, www.invitrogen.com). The integrity 
of RNA was checked on a 1.5% agarose gel 
and visualized by ethidium bromide staining. 
The concentrations and quality of RNA were 
also measured with Nanodrop ND-1000 
(Thermo Scientific, 
www.thermoscientific.com). The 260/280 nm 
absorbance ratios of all RNA samples were 
between 1.8 and 2.2. First strand cDNA was 
synthesized from 1 ug total RNA using M- 
MLV reverse transcriptase (Takara Bio Inc., 
www.takara-bio.com) and Oligo (dT18) as the 
anchor primer. The reaction mixtures were 
incubated at 70 °C for 10 min followed by 42 
°C for one hour and 70 °C for 15 min. 

Quantitative real-time PCR 

The qPCR reactions were carried out with 
SYBR Premix Ex Taq (Takara Bio Inc.) 
following the manufacturer's protocol using 
an ABI Prism 7300 (Applied Biosystems, 
www.appliedbiosystems.com). The reaction 
mixture consisted of 2 uL of cDNA template 
in a final reaction volume of 20 uL. The qPCR 
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protocol included an initial step of 95 for 30 
sec, followed by 40 cycles of 95 °C for five 
sec and then annealed at 60 °C for 31 sec, 
followed by one cycle of 95 °C for 15 sec, 60 
°C for 60 sec, and 95 °C for 15 sec. Products 
were dissolved by curve analysis (60-95 °C) 
after 40 cycles. The specificity of the qPCR 
reactions was monitored with melting curve 
analysis using SDS software (version 1 .4) and 
gel electrophoresis. Amplification efficiencies 
were determined by a series of template 
dilutions. All experiments were repeated in 
triplicate. 

Data processing 

The raw Ct values were obtained using the 
SDS software of ABI 7300 (version 1.4). The 
algorithms including geNorm (Vandesompele 
et al. 2002), Normfinder (Andersen et al. 
2004), ACt approach (Silver et al. 2006), and 
Stability index (Brunner et al. 2004) were 
used to analyze the raw Ct values of selected 
HKGs. The analysis procedures strictly 
followed the manuals of the algorithms. The 
fold changes of Siwi-\ and Siwi-2 genes in 
different developmental stages or tissues were 
calculated using standard delta-delta-Ct 
method (Pfaffi 2001). 

Results 

Identification of candidate housekeeping 
genes 

According to the mining of qPCR 
publications, six commonly-used candidate 
HKGs were selected to validate their 
expression stability, including cytoplasmic 
actin gene A3 (actin A3), cytoplasmic actin 
gene Al (actin Al), glyceraldehyde-3- 
phosphate dehydrogenase {GAPDH), 
glycerol-3 -phosphate dehydrogenase-2 
(G3PDH), E2F transcription factor 4 (E2F), 
and ribosomal protein 49 (rp49). The 
sequences of six HKGs in B. mori were 
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obtained by searching the refseq database of 
NCBI (Walter et al. 2009). The HKGs of three 
lepidopteran pests, C. suppressalis, P. 
xylostella, and S. exigua were searched from 
the transcriptome data using the BLASTX 
program (unpublished data). The sequences of 
actin A3 in B. mori were used to search for 
homologs in the three lepidopteran pests, but 
homologs were not found and thus not used in 
further experiments. Rp49 of S. exigua was 
not considered for further analysis because it 
had a very short sequence of poor quality. 

Primer specificities were confirmed by 
melting curve analysis and all primer pairs 
amplified a single PCR product with the 
expected sizes. PCR products were confirmed 
by bi-direction sequencing. Amplification 
efficiencies of primers reached the standard 
requirements of conventional qPCR, which 
were determined by a serial dilution of cDNA 
template. All the experiments including RNA 
purification, cDNA synthesis, and qPCR were 
carried out by the same individual using the 
same reagent lot. 

Expression levels of the reference gene 
candidates 

The Cycle threshold (Ct) values in qPCR gave 
an overview of the gene expression variation 
in the samples. All selected HKGs had 
moderate abundance in different 
developmental stages of C. suppressalis. The 
mean Ct values of HKGs in four insect 
species were 15.20-29.56 cycles. Individual 
HKGs had different expression levels across 
all tested samples. According to the variations 
of Ct values, GAPDH showed the smallest 
gene expression variation (below five cycles) 
in B. mori and S. exigua, while actin Al had 
the highest expression variation (above eight 
cycles). However, actin Al was the most 
stable gene C. suppressalisand P. xylostella, 
whereas GAPDH was the least stable. The 
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Figure I . The variance of raw Ct values of selected housekeeping genes (HKGs) in all tested samples. High quality figures are 
available online. 



wide range of Ct values of six tested HKGs 
indicated that there was no HKG with 
constant expression under these conditions 
(Figure 1). Thus, determining a suitable HKG 
as a reference gene requires careful 
confirmation of expression stability. 

Expression stability of selected HKGs 
across different developmental stages 

Gene expression patterns across different 
developmental stages are frequently 
investigated. Insects were collected from 
developmental stages including egg, larvae 
collected at the first day of each instar, pupa, 
and adult, in order to study the expression 
stability of selected HKGs. Raw Ct values 
were first transformed to relative quantities 
using the ACt method. Then, geNorm and 
NormFinder were used to calculate the 
expression stability of selected HKGs. 
geNorm estimates the gene expression 
stability by calculating stability measure M for 
an HKG as the average pairwise variation for 
that gene with all other tested HKGs. The 
HKGs with low M values are the stably 
expressed genes. Based on this principle, 
geNorm recommend a pair of the most stable 



HKGs as candidate reference genes 
(Vandesompele et al. 2002). According to the 
results of geNorm, the stabilities of selected 
HKGs were actin A3 = rp49 > GAPDH > 
G3PDH > E2F > actin Al in B. mori. The 
HKG stabilities in P. xylostella were Actin Al 
= E2F > G3PDH > rp49 > GAPDH; in C. 
suppressalis actin Al = G3PDH > E2F > 
rp49; in S. exigua GAPDH = G3PDH>E2F> 
actin Al. These results demonstrated that 
actin Al was most stable in P. xylostella and 
C. suppressalis, whereas it was least stable in 
the other two lepidopteran insects across 
different developmental stages. GAPDH was 
the most stable HKG in S. exigua, whereas it 
was the least stable in P. xylostella (Table 2 
and Figure SI). These results proved that the 
stabilities of HKGs are quite different in the 
four lepidopteran insects used in this study. 

To confirm the results of geNorm, 
NormFinder was used to analyze the qPCR 
data. In general, the optimum reference genes 
found by NormFinder were similar to those of 
geNorm, suggesting the reliability of data 
analysis. Both geNorm and NormFinder 
ranked actin Al as the most stable HKG in P. 
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Figure 2. ACt approach analysis of expression stability of selected housekeeping gene (HKGs) in the (A) Bombyx mori, (B) 
Plutella xylostella, (C) Chilo suppressalis, and (D) Spodoptera exigua. The qPCR data from all tested sample were put together for 
analysis. ACt variability of HKGs is shown as medians (lines), 25 th percentile to the 75 th percentile (boxes), and ranges 
(whiskers) estimated from the data of all samples. High quality figures are available online. 



xylostella but the least stable in the B. mori 
and S. exigua. The least stable HKGs found 
by both two algorithms were the same. There 
was an inconsistency in determining the most 
stable gene in B. mori and C. suppressalis by 
two algorithms. In B. mori, GAPDH was 
evaluated as the most stable gene by 
NormFinder, but rp49 and actin A3 were the 
best candidates for reference genes by 
geNorm. In C. suppressalis, E2F was the most 
stable by NormFinder whereas actin Al and 
G3PDH were the most stable HKGs (Table 
2). Taking all factors into account, we 
recommend that rp49 in the B. mori, actin Al 
in P. xylostella and C. suppressalis, and 
GAPDH in S. exigua be used as the reference 
genes across different developmental stages. 

In some cases, using the two HKGs in 
combination is more accurate than just using 
one reference gene. NormFinder not only 
finds the optimum reference genes out of a 
group of candidate genes, but also (in contrast 



to geNorm) takes information of groupings of 
samples into account (Andersen et al. 2004). 
Therefore, NormFinder is able to estimate the 
variation between sample groups, which can 
determine the best combination of two 
reference genes for normalization. The best 
combination of HKGs was GAPDH and 
G3PDH in B. mori with stability value of 
0.421. The combinations of actin Al and 
G3PDH were estimated as the best pair of 
reference genes in both P. xylostella and 5*. 
exigua. In C. suppressalis, E2F and G3PDH 
was the best combination with stability value 
of 0.454 (Table 3). It should be noticed that 
choosing the best combination considers the 
pairwise inter-variations and intra-variations 
in the grouping of samples. It does not rely on 
the stability value of an individual HKG. 
Actin Al is the least stable HKG in the S. 
exigua when analyzed singly. However, it was 
estimated to be one of the best combinations 
{actin Al and G3PDH) when considering 
multiple reference genes. 
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Supplementary Figure 2. The expression stability of selected housekeeping genes (HKGs) in different tissues estimated by 
geNorm. High quality figures are available online. 



Expression stability of candidate reference 
genes in different tissues 

The gene expression pattern in different 
tissues is also frequently studied to infer the 



gene functions. To screen suitable reference 
genes, different tissues were collected 
including head, midgut, ovary, testis, fat body, 
malpighian tube, epidermis, or silk gland of B. 
mori, C. suppressalis, and S. exigua. The 
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expression stabilities of selected HKGs were 
analyzed with geNorm and NormFinder. 
Though there were differences between the 
ranks of HKGs by two algorithms, the most 
stable and the least stable HKGs in different 
tissues of three insects found by NormFinder 
were identical with those found by geNorm. 
Considering all factors, GAPDH was the most 
stable of HKGs in B. mori, while actin Al in 
C. suppressalis and E2F in S. exigua were the 
best candidate reference genes (Table 4 and 
Figure S2). In C. suppressalis, the expression 
of actin Al was stable among different 
developmental stages and tissues. However, 
the most stable HKGs in the different tissues 
of B. mori and S. exigua were not same as 
those across different developmental stages, 
suggesting that HKGs have varied expression 
levels under different conditions. 

For the best combinations of two HKGs, 
GAPDH and G3PDH were the most suitable 
candidates in different tissues of B. mori, 
which was identical with those across 
different developmental stages. Rp49 and 
actin Al in C. suppressalis and GAPDH and 
E2F in S. exigua were the best combinations 
in different tissues, which was different from 
those across different developmental stages 
(Table 3). 

Expression stability of candidate reference 
genes in all tested samples tissues 

It is impossible to evaluate the expressions of 
HKGs under all kinds of experimental 
situations. To provide a general understanding 
of the expression stability of selected HKGs, 
the qPCR data from all tested samples was put 
together for evaluation. According to the 
results of geNorm and NormFinder, the most 
stable gene in B. mori was GAPDH followed 
by actin A3. Actin Al was the least stable 
HKG in B. mori and S. exigua, but it was the 
most suitable candidate of reference genes in 
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C. suppressalis and P. xylostella. GAPDH was 
the most stable HKG in S. exigua, followed by 
G3PDH (Table 5). 

To validate the analysis results of geNorm and 
NormFinder, stability index assay and ACt 
analysis were also used to estimate the 
expression variations of HKGs in all tested 
samples. The stability index assay was first 
used to select suitable internal controls during 
the development of poplar (Brunner et al. 
2004). According to the stability index assay, 
rp49 had the lowest stability index in B. mori, 
P. xylostella, and C. suppressalis, while 
GAPDH was the most stable HKG in S. 
exigua (Table 5 and Table SI). In general, 
stability index assay had different results with 
those of geNorm and NormFinder in 
identifying the most stable HKGs. However, 
the most stable HKGs found by geNorm and 
NormFinder also had a low stability index. 

The ACt approach was employed to select the 
most suitable HKG in reticulocytes (Silver et 
al. 2006). The principle of the ACt approach 
was to examine the ACt value between the 
two HKGs in different samples. If the ACt 
value remains constant, it suggests that either 
both HKGs are stably expressed among the 
tested samples or are co-regulated. However, 
the fluctuation of the ACt value indicates that 
at least one of them is variably expressed. 
After introducing more HKGs into the 
comparisons, an appropriate HKG can be 
found for a particular experimental system 
(Figure 2). The most stable HKGs selected by 
the ACt approach were identical with those by 
geNorm and NormFinder, which were 
GAPDH in B. mori and S. exigua and actin Al 
in P. xylostella and C. suppressalis (Table 5 
and Table S2). 
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The impact of reference genes on qPCR 
data analysis 

Based on the above analysis, we present 
evidence that the expression of HKGs in 
lepidopteran insects was not constant as 
expected. A stable HKG in one species did not 
mean that it also was expressed stably in 
another. Thus, it is our interest to study the 
impact of using "wrong" reference genes on 
the qPCR data analysis. The Siwi-l and Siwi-2 
genes in B. mori were chosen for evaluations 
of different HKGs. The samples from 
different developmental stages and six 
different tissues of the fifth instar larvae were 
used in the experiments. 

For qPCR data analysis across different 
developmental stages, the egg stage was 
designated as the "calibrator sample" 
(abundance set to IX). The relative expression 
levels of all other samples were then 
calculated relative to the calibrator sample. 
The results indicated that using different 
reference genes had apparent effects on the 
qPCR analysis. Rp49 was recommended to be 
the most suitable candidate of reference genes 
across different developmental stages. The 
relative abundance of both Siwi-\ and Siwi-2 
using rp49 as the reference gene were in 
general similar to those obtained when using 
actin A3, the second most stable HKG (Table 
6). However, the results normalized using 
rp49 were significantly different from those 
using actin Al, the least stable HKG, as the 
internal control. In the third instar larvae, the 
relative abundance of Siwi-l was 0.314 when 
using rp49 as the reference gene, but it 
changed to 0.067 when using actin Al. A big 
difference could also be observed for the Siwi- 
2 gene. The relative abundance of Siwi-2 
normalized using actin Al (about 10~ 4 level) 
was significantly lower than those using rp49 
(about 10" 2 level) (Table 6). 



Teng et al. 

For qPCR data analysis in different tissues, 
the head was selected as the "calibrator 
sample". In this case, GAPDH was the most 
stable HKG and was suggested to be the best 
reference gene. If using GAPDH as the 
reference gene, the relative mRNA 
abundances of Siwi-l in the ovary and testis 
were only 2.406 and 1.769, respectively, but 
they changed to 26.438 and 46.709 when 
using actin Al as the internal control, which 
was the least stable HKG (Table 7). For Siwi- 
2 gene, if GAPDH was used as the reference 
gene, the relative expression levels in ovary 
and testis were 1.214 and 2.488, respectively. 
But they were 13.347 and 64.792 when using 
actin Al for normalization. On the contrary, 
the results of GAPDH were similar with those 
of actin A3, the second stable HKG (Table 7). 

These surprising differences demonstrate the 
importance of choosing a proper reference 
gene. If an unsuitable reference gene were 
used in normalizing qPCR data, incorrect 
results might be obtained and lead to improper 
conclusions. 

Discussion 

qPCR has become a powerful tool for 
estimating gene expression levels because of 
its sensitivity and accuracy (Bustin 2000; 
Gachon et al. 2004; Wong and Medrano 2005; 
Nolan et al. 2006). The variations of qPCR 
will be unavoidably introduced by the 
procedures of RNA preparation, cDNA 
synthesis, or PCR processing. These 
variations can be controlled to some extent by 
carefully conducting the qRT-CPR 
experiments. The Ct values used for qPCR 
data analysis is determined mostly by the 
quantity of cDNA template (Walker 2002). 
The cDNA quantity is determined by the total 
RNA used in the first strand cDNA synthesis 
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and enzymatic efficiency of reverse 
transcriptase. Therefore, there is a close 
relationship between RNA integrity and 
expression quantitation. The quality and 
purity of the total RNA should be stringently 
analyzed before use in qPCR analysis (Bustin 
and Nolan 2004; Fleige and Pfaffl 2006; 
Fleige et al. 2006). However, it is impossible 
to eliminate all variations by these strategies. 
It is necessary to use HKGs as the proper 
reference genes for normalization. 

To be a good reference gene, an HKG should 
meet three criteria: first, it should have 
amplification efficiency similar to the target 
genes; second, it should have moderate 
expression level; and third, its expression 
should be stable in all test samples. 
Unfortunately, it is unlikely that a single 
universal reference gene can be found in all 
species or in all experimental conditions of a 
species (Bustin 2000). Almost all genes 
including HKGs are regulated by other 
"regulators". It is unsafe to deduce the 
expression stability of an HKG based on 
known information of gene regulation. The 
actin gene is a widely used reference gene 
because it is thought to be stably expressed in 
all possible conditions. However, this gene 
has been challenged for its suitability as the 
internal control (Selvey et al. 2001). In our 
work, though actin Al was the most stable in 
P. xylostella and C. suppressalis, its 
expression was the least stable in B. mori and 
S. exigua. 

Evaluation of the expression stability requires 
mathematical methods. Many algorithms such 
as geNorm (Vandesompele et al. 2002), 
Normfinder (Andersen et al. 2004), ACt 
approach (Silver et al. 2006), and stability 
index (Brunner et al. 2004) have been 
developed. The most stable HKG found by 
geNorm, NormFinder, and ACt approach were 
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the same in all four lepidopteran insects. But 
stability index had different results. We 
reasoned that this approach only considers the 
stability of HKGs individually but the other 
three algorithms evaluate the stability of 
HKGs by analyzing its variance comparing 
with other HKGs in pair. Thus, the reliability 
of stability index is not as high as others. 
Although the expression stability analyzed by 
ACt approach was similar to that by geNorm 
and Normfinder, it is not a good choice for 
reference gene selection because it relies on a 
strict principle that the expression variance of 
two HKGs should be identical in all samples 
in all experimental conditions or cell types. In 
addition, the computation procedures of ACt 
approach are cumbersome (Silver et al. 2006). 
geNorm and NormFinder use different 
mathematical methods to estimate the 
expression stability. The results of two 
algorithms can be used for cross validation. 
geNorm is one of the reliable algorithms to 
search for stably expressed HKGs that have 
low intra group variation and non-vanishing 
inter-group variation (Jain et al. 2006; 
Exposito-Rodriguez et al. 2008; Gutierrez et 
al. 2008; Jian et al. 2008). Compared with the 
geNorm, NormFinder finds not only the most 
stable gene but also the best combination of 
two reference genes if multiple internal 
controls are considered in qPCR data analysis 
(Oelkers et al. 2008). 

In summary, our work presented a number of 
the most stable HKGs that are suitable to be 
used as the reference genes in four 
lepidopteran insects. The results demonstrated 
that no single universal reference gene could 
be used in all experiment conditions or 
treatments. The most reliable approach of 
choosing suitable HKGs as reference genes is 
to validate their expression stability using 
algorithms such as geNorm and NormFinder. 
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Table I . The primer sequences of selected housekeeping genes (HKGs) used in the experiments. 



Species 


Genes 


Cenbank Accession No. 


Primers (from S'to 3') 


Sense 


Antisense 


Bombyx mori 


Actin A3 


X04507.1 


TACGCACTCCCTCACGCCATCC 


CCTCTGGGCAACGGAATCTTTCG 


actin Al 


NM 001126252.1 


GAGAGGGAAATCGTTCGTGACATC 


AGGAAGGAAGGCTGGAAGAGAG 


GAPDH 


AB262581.1 


CGCTGGAATTTCTTTGAATGAC 


CAATGACTCTGCTGGAATAACC 


G3PDH 


AD164061.1 


GCCTTCTAATGTTGTTGCTGTTCC 


ACCACCTTCGGCTATATCAAATGC 


E2F 


DQ311161.1 


ATTCAATGGAGAGGAGCAGGTC 


CATATCTTGGTTGTCTGGTTCGTC 


rp49 


AD048205.1 


AGGCATCAATCGGATCGCTATG 


TTGTGAACTAGGACCTTACGGAATC 


Plulella xytostella 


actin Al 


Transcriptomc data 


CACGGCA1 1 G 1 CACAGA1 1 GGAAC 


CACTTCGGCGGAC Tl CTCACG 


GAPDll 


CCGCATGGTATGTTCATTGAC 


CAGGACTGTTCACTACTCTCG 


GiPDH 


CCTACTGCTGCTGCGTTGTC 


GCCACCTCAGAAGCGATGTTG 


E2F 


TCTCAACTTCCTCGTCTTCCATAG 


ACCAACAACCAAGTATTAGTGCTG 


ip49 


CCTGGTTCCTGACGCTGTTC 


GCATCGCCTGTCCATCTTCC 


Chilo suppressalis 


actin Al 


Transcriptome data 


GTCGCTTCCCAAATTACATC 


CTCCATATCGTTCCAGTCG 


GAPDH 


ATCTTGTGACCATCAATG 


GCTGAATATGCTGCTTAC 


G3PDH 


GTTGTGCCTCACCAATTTGTCAG 


GCCACCTTCAGCGATGTCG 


E2F 


ATTGCTGTGTGATAAAGAAGAAC 


AGAAGGTGGTGGACTCAAC 


ip49 


TGTTAGACACCATTCAGATAGG 


GACCACAAGGAAGAAGATGC 


Spodoptera exigua 


actin Al 


Transcriptomc data 


TTTCTCACGGTTCGCTTTGG 


GTATCCTCACGCTCAAGTATCC 


GAPDH 


AACATTTATCTCTACAACGCAATC 


GTGACAACCACTCATCTATCTTC 


GiPDH 


GTAGTTTCCCACCAGTTTGTCAG 


CACCTCAGACGCAATGTTAGC 


E2F 


GCCCAATCTGTTTCCTCAAAAG 


GAACTTGCTCGCCGTAAGAC 



Table 2. The stability values of selected housekeeping genes 
(HKGs) across the different developmental stages. 







Norm Finder 


geNorm 


Species 


Genes 


Stability 
values 


Rank 


Stability 
values 


Rank 




Actin A3 


0.579 


3 


0.385 


1 




actin Al 


1.085 


6 


1.547 


6 


Bombyx mori 


GAPDH 


0.459 


1 


0.548 


3 


G.IPDH 


0.743 


4 


0.897 


4 




E2F 


0.752 


5 


1.27 


5 




rp49 


0.52 


2 


0.385 


1 




actin A I 


0.428 


1 


0.417 


1 




GAPDH 


1.081 


5 


1.532 


5 


Plulella xyloslella 


GiPDH 


0.685 


2 


0.968 


3 




E2F 


0.751 


3 


0.417 


1 




rp49 


0.781 


4 


1.311 


4 




actin A 1 


0.647 


2 


1.087 


1 


Chilo suppressalis 


G3PDH 


0.663 


3 


1.087 


1 


E2F 


0.505 


1 


1.192 


3 




rp49 


0.917 


4 


1.371 


4 




actin Al 


2.012 


4 


3.36 


4 


Spodoptera exigua 


GAPDH 


1.024 


1 


1.869 


1 


G3PDH 


1.19 


2 


1.869 


1 




E2F 


1.454 


3 


2.287 


3 



Table 3. The best combination of two housekeeping genes (HKGs) recommend by NormFinder, which can be used as the 
multiple references genes in different developmental stages or tissues. 



Species 


Developmental stages 


Tissues 


Best combination 


Stability value 


Best combination 


Stability value 


Bombyx mori 


GAPDH and GiPDH 


0.421 


GAPDH and G3PD1I 


0.532 


Plutella xyloslella 


actin A 1 and GiPDH 


0.474 


N/A 


N/A 


Chilo suppressalis 


E2F and GiPDH 


0.454 


rp49 and actin Al 


0.722 


Spodoptera exigua 


actin Al and GiPDH 


1.185 


GAPDH and E2F 


0.621 
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Table 4. The stability values of selected housekeeping genes 
(HKGs) in different tissues. 







NormFinder 


geNorm 


Species 


Gene name 


Stability 
values 


Rank 


Stability 
values 


Rank 




act in A3 


0.875 


3 


0.568 


1 




actinAl 


1.847 


6 


2.03 


6 


Bomhyx tnori 


GAPDH 


0.364 


1 


0.568 


1 


G3PDH 


1.011 


5 


1.29 


5 




E2F 


0.923 


4 


1.123 


4 




rp49 


0.822 


2 


0.646 


3 




actinAl 


0.475 


1 


1.114 


1 




GAPDH 


1.503 


5 


2.774 


5 


Chilo suppressalis 


G3PDH 


1.061 


2 


1.114 


1 




E2F 


1.315 


4 


1.653 


4 




rp49 


1.136 


3 


1.3 


3 




actinAl 


1.419 


4 


1.732 


4 


Spodoptera exigua 


GAPDH 


0.673 


1 


1.182 


3 


G3PDH 


0.831 


3 


0.598 


1 




E2F 


0.673 


1 


0.598 


1 



Table 5. The stability values of selected housekeeping genes (HKGs) in all tested 
samples. 







Norm finder 


GeNorm 


Stability 


index 


ACt 


Species 


Genes 


Stability 
values 


Rank 


Stability 
values 


Rank 


Stability 
values 


Rank 


SD 


Rank 




Actin A3 


0.646 


2 


0.595 


1 


0.364 


2 


2.202 


3 




act in Al 


1.587 


6 


1.89 


6 


3.451 


6 


3.439 


6 


Bombyx tnori 


GAPDH 


0.555 


1 


0.595 


1 


0.409 


3 


1.856 


1 


G3PDH 


0.91 


4 


1.56 


5 


1.803 


5 


2.325 


5 




E2F 


0.854 


3 


1.369 


4 


1.562 


4 


2.281 


4 




rp49 


1.042 


5 


0.877 


3 


0.083 


1 


2.095 


2 




actin A 1 


0.428 


1 


0.649 


1 


3.565 


2 


1.202 


1 




GAPDH 


1.081 


5 


1.588 


5 


3,958 


4 


2.235 


5 


Plutella xylostella 


G3PDH 


0.685 


2 


0.923 


3 


3.704 


3 


1.498 


3 




E2F 


0.751 


3 


0.649 


1 


5.585 


5 


1.439 


2 




rp49 


0.781 


4 


1.156 


4 


2.232 


1 


1.6 


4 




actin A I 


0.479 


1 


1.47 


1 


1.404 


2 


1.84 


1 




GAPDH 


1.838 


5 


2.339 


5 


11.79 


5 


3.378 


5 


Chilo suppressalis 


G3PDH 


1.269 


4 


1.936 


4 


4.706 


4 


2.625 


4 




E2F 


0.846 


2 


1.47 


1 


1.799 


3 


2.361 


3 




rp49 


1.027 


3 


1.766 


3 


1.338 


1 


2.197 


2 




actin Al 


2.426 


4 


2.844 


4 


6.454 


4 


2.592 


4 


Spodoptera exigua 


GAPDH 


0.443 


1 


1.765 


1 


2.249 


1 


1.8 


1 


G3PDH 


1.058 


2 


1.765 


1 


3.962 


2 


1.913 


2 




E2F 


1.187 


3 


1.939 


3 


3.176 


3 


2.264 


3 



/ 

Table 6. The relative mRNA abundance of Siwi-\ and S/vw-2 genes across different 
developmental stages of the silkworm when using different housekeeping genes (HKGs) 
as the internal control. The rank of reference genes was determined by geNorm. 



Genes 


Rank 


r P 49 


actin A3 


GAPDH 


G3PDH 


E2F 


actinAl 


1 


2 


3 


4 


5 


6 


Siwi-l 


Egg 


1 


1 


1 


1 


1 


1 


1" 


0.169 


0.272 


0.202 


0.316 


0.952 


0.037 


2 nt * 


0.087 


0.12 


0.146 


0.195 


0.278 


0.016 


3 rd 


0.314 


0.327 


0.292 


0.084 


0.407 


0.067 


4* 


0.027 


0.025 


0.041 


0.128 


0.155 


0.293 


5* 


0.048 


0.063 


0.083 


0.088 


0.179 


0.01 


Pupa 


0.773 


1.55 


0.498 


0.379 


0.604 


0.157 


Adult 


0.855 


0.942 


0.625 


0.693 


0.073 


0.245 


Siwi-2 


Egg 


1 


1 


1 


1 


1 


1 


r 


0.004 


0.031 


0.022 


0.175 


0.577 


0.00001 




0.173 


0.166 


0.879 


0.02 


2.703 


0.00029 




0.019 


0.021 


0.066 


0.011 


0.479 


0.00003 


4* 


0.033 


0.039 


0.178 


0.015 


1.162 


0.00005 


5* 


0.028 


0.032 


0.087 


0.031 


0.937 


0.00003 


Pupa 


0.194 


0.292 


0.983 


0.147 


1.994 


0.00044 


Adult 


0.209 


0.109 


0.366 


0.058 


1.198 


0.00015 
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Table 7. The relative mRNA abundance of S/vw-l and Siwi-2 genes in different tissues of 
the silkworm when using different housekeeping genes (HKGs) as the internal control. 
The rank of reference genes was determined by geNorm. 





IV a LI K 


GAPDH 


actin A3 


rp49 


E2F 


(13PDH 


actin Al 


1 


2 


3 


4 


5 


6 




Head 


1 


1 


1 


1 


1 


1 




Midgut 


0.498 


0.151 


5.312 


0.815 


0.218 


1.962 


Siwi-l 


Ovary 


2.406 


1.599 


2.037 


0.532 


2.513 


26.438 


Testis 


1.769 


1.41 


1.153 


0.255 


0.701 


46.079 




Silk gland 


0.788 


0.176 


0.256 


1.977 


5.955 


11.864 




Malpighian tube 


1.215 


0.588 


1.474 


1.427 


0.328 


13.664 




Head 


1 


1 


1 


1 


1 


1 




Midgut 


0.827 


0.25 


8.816 


1.352 


0.362 


3.257 


Siwi-2 


Ovar>' 


1.214 


0.807 


1.028 


0.269 


1.269 


13.347 


Testis 


2.488 


1.983 


1.621 


0.358 


0.985 


64.792 




Silk gland 


0.625 


0.14 


0.203 


1.569 


4.724 


9.411 




Malpighian tube 


1.415 


0.685 


1.718 


1.663 


0.383 


15.923 



Supplementary Table 1 . Evaluation of the stability of selected housekeeping genes (HKGs) in all tested 
samples using stability index assay. HKGs with the lowest stability index are treated as the best control. 




S|M.*Cic* 


Gene names 


Mean 


Sldev 


Slope 


R 1 


CV°/. A 


Stability inde\ " 


rank 


Bombyx 
morl 


actin A3 


15.337 


0 672 


1 1 1 S • 


0 961 


4 3t 


1 3, I 


s 


actin Al 


15 634 




0 248 


0812 


13 V45 


3 451 


6 


i, vi on 


16 


;i 7; S 


i IIX'i 


1 1 Ji 1 7 


I I.I 3 


1 I'M 


3 




cr-pnii 


19884 


1 7,„, 


021 


1 'I5„ 


S 57'l 


1 S"3 


5 






21.208 


| 7,,<l 


i I'M 


0 SI 


K 06 


1 567 


4 




■ si'. 


14 X5X 


0 315 


o 039 


0 972 


2 122 


0 0X3 


1 


Pfataflts 
xytosteth 


actin A 1 


21.893 


1 840 


o 422 


0 913 


X 144 


3 565 


2 


GAPDH 


29 561 


2 15? 


0 544 


0935 


7 2X2 


3 958 


4 




G3PDH 


21 747 


1 86 


0 433 


0 949 


X 553 


3 704 


3 




t 




2 3X1 


0 54X 


0 925 


10 19J 


5 5X5 


5 




tp49 


26 075 


1.611 


0 36 


II S(,4 


I, \>7 




1 


Chtlo 

suppressahs 


actin Al 


27X54 


1 671 


0 234 


0 9X 


! 1 


1 404 




UAPDH 


?U Vs| 


1 917 


1 1165 


l 7',., 


1 1 '17*. 


1 1 7o 






cr-pnii 


21 795 




l l>5 


I 'PI 


10 3 IV 


I 7,,,, 


1 




F 1 F 


V 1 p 


1 7< t 


1 ^7 


l S3s 


7'P3 


1 -','.) 


3 




rp4V 


r, 045 


1 548 


0 199 


0 837 


6 719 


1 33X 


1 


Sptxhpfera 
rxtgiia 


actin Al 


23 675 


3 422 


0 447 


0 851 


14 455 


6 454 


1 


GAPDH 


1.: II'. 


1 705 


0 213 


0 777 


Id J7J 


2 249 




G3PDH 


20 772 


2 44 


0 337 


0 955 


1 1 746 


3 962 




E2F 


23 523 


2 348 


1 1 3 1 X 


0 919 


9 9X1 


3 |76 
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Supplementary Table 2. The raw data of gene expression stability analysis by ACt 
approach. 



Speck** 




■ 111 All u V 1 


SO 


Mm 
nimi ■ if 




ac4mA3 vs actin Al 


4.85 


3.658 


2 202 




actmA3 vs OAPDH 


■0.201 


1 285 






1 1 1 : >H 


-4.918 


2 2X6 






actaA3 vs E2F 


■6.259 


2 197 






ectaA3 vs rp49 


0 X99 


1 587 






■clin A 1 vs .x'tinA * 


0.85 


3 658 


3.439 




acta A 1 vs OAPDH 


0.649 


3.314 






■otinAl v S 03PDH 


•4.068 


^ ut'i'j 






acta Al vs E2F 


•5.409 


3 761 






actin A 1 vs rp49 


1.749 


3 395 






OAPDH vs octtnA3 


0 2*1 1 


1 285 


1 K56 




OAPDH vs octaAl 


■0.649 


3 314 






OAPDH vs 03PDH 


-4.717 






Bombyx mori 


OAPDH vs E2F 


4.058 


1 66 




OAPDH vsrp49 


1.1 


0.983 






03PDH vs uctinA.1 


4 918 


2 2X6 


2 325 




03PDH vs actinAl 


■i II6X 


i 1 169 






03PDH vs OAPDH 


4 717 


2 035 






(,iim)H .. . ;i 


•1.341 


1.756 






03PDH vsip49 


5 817 


2 477 






E2F vs actin A3 


6 259 


2 197 


2 281 




E2F vs iKlm A 1 


S 41N 


3 761 






E2F vs OAPDH 


6.058 


1.66 






1:21 5 i ' -l ! HI 


1 14 1 


1 756 






E2F vs rp49 


7 158 


2 029 






ip49 vs actin A3 


4.899 


1 587 


2 095 




rp49 vs actin A 1 


•1.749 


3 395 






rp49 vs OAPDH 


-1.1 


0 9X5 






OAPDH vs rp49 


3565 


2 215 


2 235 




OAPDH vs 03PDH 


7 929 


2.27 






OAPDH vs E2F 


6.45 


2 329 






OAPDH vs ociinAl 


7 X49 


2 126 






ip49vs OAPDH 


-3.565 


2 215 


1 6 




rp49 vs 03PDH 




1 561 








2 723 


1 522 






rp49 vs actin A 1 


4 1X2 


1 102 






03PDH vs OAPDH 


•7.929 


2.27 


1 498 


t't„r -It., 
r (Hh lid 


03PDH vs rp49 


-4.328 


1 561 




xyiosleiia 


03PDH vs E2F 


•1.605 


1.244 






03PDH vsoctinAl 


•0.146 


0.916 






E2F vs OAPDH 


■6.45 


2 329 


1 439 




E2F vs rp49 


•2.723 


1 522 






E2F vs 03PDH 


1 605 


1 244 






E2F vs acta Al 


1 IV! 


n,.f,; 






acta A 1 vs OAPDH 




2 126 


1 MS 




actin A 1 vi ip49 


-4 MB 


I 103 






actinAl vs 03PDH 


i i-. 


0.916 






actinAl vsE2F 


•1.459 
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c \ 





OAPDH vsip49 


3 52 


) 101 


3 37X 




OAFDB vs E2F 


i 5H7 


i 624 






OAPDH vs actmAl 


-2.31 


2 766 






OAPDH vsOJPDH 


? 4X6 


4 02 






rp49 vs OAPDH 


-3 52 


3 101 


2 197 




rp49vs(j3PDH 


0 123 


2 167 






rp49 vs E2F 


0 57 


2 (U4 






ip49 vs actiii A 1 


-5 29X 


1 5(15 






<i3l>nil v^CAI'IMI 




- ' 


2 625 


Cfhlo 


(,l|M)ll ■. • r 7 4 . 


-0.123 


2 167 




suppressaiis 


03PDHvsK2F 


n ' 


2 515 






03PDH vsnctinAl 


5421 


1 X 






i :i •,- «..M ! i)M 


•1 507 


3 624 


2 361 




E2Fvsrp49 


-0.57 


: ni4 






1:21 VMnl' Ml 


•0 447 


: 515 






E2F vs Dclm Al 


-5 868 


1 2X9 






actmAl vs OAPDH 


2 31 


: 766 


1 X4 




actm Al vs rp49 


5 29X 


1 5H5 






actmAl vs 0311)11 


5 421 


x 






ac4in Al vs E2F 


i X'.X 


1 'X'. 






(.API Ml v. F.2F 


-7429 


: on 


- 




OAPDH vs03PDH 


1 1 • 


1 45 






OAPDH vs actmAl 


•8 417 


1 94 






03PDHvs OAPDH 


4 H4 


1 45 


1 " ; 




03PDHvs E2F 


: 


1 62 1 




Spodoptera 


dUMMI ■ .,.ir \I 


i m 


? M.'l 




extgua 


E2F vs (JAPDH 


7 429 


: nil 


: 164 




1 :\ :-,^M< Ml 


3 3X9 


1 621 






E2F vs actm Al 




3 16 






actmAl vs OAPDH 


K 4]7 


1 94 


2 542 




actmAl vsO.IPDH 


4 569 


: 676 






actmAl vs E2F 


••NX 


3 16 
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